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ABSTRACT 

A  Loosely  Couple  Blade  Row  (LCBR)  numerical 
method  is  developed  to  analyze  the  unsteady  flow 
field  of  a  three-dimensional  multi -blade-row 
turbomachinery  problem.  This  method  allows  the  use 
of  initial  blade  numbers  for  unsteady  calculation  and 
only  one  blade  channel  per  blade  row  is  used. 
Circumferential  average  approach  to  obtain  a 
converged  steady  solution  as  the  initial  condition  is 
adopted  for  unsteady  calculation.  Relatively,  it  not 
only  maintains  the  blade  configuration,  but  is  also 
computationally  more  efficient.  The  numerical 
method  used  is  a  compressible  viscous  finite  volume 
algorithm  solving  Reynolds  averaged  Navier-Stokes 
equations.  Artificial  dissipation  terms  similar  to  that 
used  by  Jameson  is  adopted  to  suppress  the  numerical 
oscillations.  Four-Stage  Runge-Kutta  scheme  is  used 
to  advance  the  flow  equations  in  time.  Residual 
smoothing  and  multi -grid  techniques  are  employed  to 
accelerate  the  convergence.  Baldwin  and  Lomax 
algebraic  turbulent  model  is  used  for  the  calculation 
of  eddy  viscosity.  The  UTRC  (United  Technologies 
Research  Center)  large  scale  turbine  is  used  here  as 
the  test  case.  The  unsteady  flow  predicted  correlates 
well  with  data.  The  results  indicate  that  the  unsteady 
turbomachinery  flowfield  can  be  solved  with  single 
flow  channel  of  unequal  pitch  on  each  blade  row  . 

LJNTRQDUCUQJ^ 

Multistage  turbomachinery  calculations  have 
becoming  popular  since  mid  80’s.  The  two  main 
streams  are  circumferential  average  calculations  and 
unsteady  calculations.  Among  the  circumferential 
average  methods,  representatives  studies  are  from  Ni 
[1989],  Adamczyk  [1990],  Denton  [1992]  and  Dawes 
[1992].  These  methods  have  been  used  to  replace  the 
traditional  row-by-row  design  methods  used  in 
industry  for  many  years.  However,  these  type  of 


methods  generally  ignore  the  unsteady  interactions 
between  blade  rows.  Despite  this  kind  of  method 
ignores  true  unsteady  interaction  effects  between 
blade  rows,  it  allows  the  engineers  to  design  and 
optimize  a  multi-stage  turbomachine  within  a  much 
shorter  time.  The  advantage  of  running 
circumferential  average  methods  is  more  obvious  as 
the  number  of  blade  rows  is  increased. 

The  unsteady  pressure  fluctuations  on  blade 
surfaces  can  be  significant  especially  for  modem  high 
loading  blade  designs  with  smaller  gaps  between 
blade  rows.  As  a  result,  other  means  to  calculate  the 
fluctuating  pressures  on  the  blades  must  be  found  in 
order  to  assure  the  design  fatigue  life.  An  approach 
using  three-dimensional  quasi-steady  methodology 
was  studied  by  Chen  and  Lee  [1998].  It  does  provide 
pressure  fluctuation  data  obtained  from  this  method 
and  is  more  econmomical  than  full  unsteady  analysis. 
However  this  can  only  be  considered  as  an 
intermediate  approach.  The  simulation  of  fully 
coupled  unsteady  rotor-stator  interactions  allows  the 
total  understanding  of  the  flow  physics.  The  three- 
dimensional  fully  coupled  unsteady  rotor-stator 
interaction  simulation  starts  from  Rai  [1987]  with  an 
overlaid  patched  grid  system.  Other  methods 
thereafter  (e.g.  Hah  [1992])  more  or  less  use  the 
similar  approach.  Their  calculations  frequently 
modified  the  rotor  and  stator  numbers  to  a  simple 
arithmetic  ratio  like  2:3  or  1:1,  The  change  of  blade 
number  was  accompanied  by  a  change  of  blade  size 
in  order  to  nnaintain  its  solidity  or  blockage.  This 
however  biased  the  physical  problem.  The  predicted 
amplitude  and  frequency  in  the  rotor-stator  interaction 
are  also  biased.  Moreover,  these  methods  are  very 
computational  costly.  As  a  result,  they  were  mostly 
used  in  the  final  design  verification  or  when  a 
problem  occurs.  Erdos  et  al.  [1977]  proposed  a  phase- 
lagged  method  for  unsteady  analysis.  This  method 
uses  only  one  blade  pitch  for  each  blade  row.  It 
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requires  to  store  the  unsteady  flow  data  over  an  entire 
period  so  that  the  blades  can  feel  periodical 
excitations  from  neighboring  blade  rows.  The 
memory  requirement  is  thus  tremendous.  This 
situation  will  become  worse  as  the  number  of  blade 
rows  is  large. 

A  simplified  loosely  coupled  approach  for 
unsteady  turbomachinery  flowfield  predictions  was 
investigated  by  Hodson  [1985]  and  Ho  and 
Lakshminarayana  [1998].  They  simply  model  the 
upstream  unsteady  disturbances  with  an  incoming 
wakes.  The  wakes  can  be  calculated  separately  from  a 
steady  state  solution  from  an  upstream  blade  row. 
Domey  et  al.  [1996]  developed  an  iterated  loosely 
coupled  blade  row  method  and  applied  it  to  a  two- 
dimensional  turbine  stage.  In  their  approach,  the 
perturbations  between  two  adjacent  blade  rows  were 
iterated  to  updated  values  until  convergence.  It  shows 
that  the  loosely  coupled  calculation  can  predict 
unsteadiness  quite  well  but  is  an  order  of  magnitude 
in  computational  efficiency  compared  with  a  fully 
couple  unsteady  analysis  in  a  two-dimensional  sense. 
Ho  and  Lakshminarayana’s  three-dimensional  loosely 
coupled  analysis  modified  the  rotor/stator  blade 
numbers  from  42:40  to  40:40  so  an  equal  pitch 
channel  for  stationary  and  rotating  blade  rows  can  be 
used.  Blade  size  were  scaled  to  maintain  its  blockage 
effect.  Domey  et  al.  modified  the  stator/rotor  numbers 
from  22:28  to  21:28  so  an  3:4  blade  ratios  can  be  used 
to  have  a  periodic  boundary  condition  satisfied. 
Neither  of  the  above  loosely  coupled  method  avoid 
the  blade  geometry  modification. 

The  purpose  of  the  present  study  is  to  develop  a 
method  that  can  solve  a  turbomachinery  unsteady 
flowfield,  but  also  has  the  following  distinct 
advantages:  Use  only  one  blade  channel  for 
calculation  while  keeping  the  blade  configuration 
unchanged.  The  numerical  method  used  in  the  present 
study  is  a  compressible  viscous  finite  volume 
algorithm  solving  Reynolds  averaged  Navier-Stokes 
equations.  The  artificial  dissipation  terms  similar  to 
that  used  by  Jameson  is  adopted  to  suppress  the 
numerical  oscillations.  Four-stage  Runge-Kutta 
scheme  is  used  to  advance  the  flow  equations  in  time. 
Residual  smoothing  technique  is  employed  to  further 
accelerate  convergence.  Baldwin  and  Lomax 
algebraic  turbulent  model  [1978]  is  used  for  the 
calculation  of  eddy  viscosity.  The  UTRC  large  scale 
turbine  test  data  were  used  to  validate  the  CFD 
calculations. 


2.  NUMERICAL  METHOD 

The  governing  three-dimensional,  compressible 
Navier-Stoke  equations  in  integral  form  can  be 
written  as: 

—\WdV+^Fc-dA=^F^-dA+\HdV  (1) 

dt 

where 

F^-El-\-Fj-^Gk 
Fy  =Ri  +S  J'^T  k 
and 
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where  e  is  the  total  energy  per  unit  volume, 
v'=  {u\v\w')  is  the  relative  velocity,  its  relationship 

with  the  absolute  velocity  V  is: 

v'  =  v  -i2xF 

,  F  =  (x,y,z) 

The  algebraic  two-layer,  eddy  viscosity  turbulence 
model  derived  by  Baldwin  and  Lomax  is  used  for  the 
calculation  of  turbulent  eddy  viscosity. 

The  physical  domain  is  divided  into  many  small 
control  volumes.  The  discretization  form  of  eq.  (1)  is 

d 

—{hw)+Q^-^Qy=hH  (2) 

dt 

where  h  is  the  volume  of  grid,  and 
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Qc=^^ck'^k 

k 

Qv  =^^vk  -^k 
k 

where  k  =1-6,  S  is  the  surface  area  on  each  side  of 
control  volume. 

Equation  (2)  is  solved  using  central  differencing 
scheme  which  is  second  order  accurate  in  space.  In 
order  to  damp  the  numerical  oscillation  associated 
with  central  differencing  scheme,  artificial  dissipation 
terms  similar  to  that  used  by  Jameson  [1981]  are 
considered.  Eq.  (2)  can  be  written  as: 

—{hw)+Q^^Q^-D=hH  (3) 

dt 

where  D  is  the  dissipative  term: 

D=(D^  +£)^  )w  (4) 


wherea=2/3n  is  the  time  step  in  the  ^-direction. 
It  can  be  approximated  as: 


Four-stage  Runge-Kutta  scheme  is  used  to  advance 
the  flow  equation  (3)  in  time.  To  further  accelerate 
convergence  and  allow  larger  time  steps  be  used, 
implicit  residual  smoothing  technique  is  used  to 
minimize  numerical  instability.  The  residue  R 
calculated  above  is  smoothed  after  each  Runge-Kutta 
stage  by: 

(10) 

where  standard  second  order 

differencing  operator,  and  smoothing  parameters 
s ^  =£^  =  £^=0.5  are  chosen  . 


The  ^“direction  operator  is  given  by: 

w-C^  ( ^2  -F4  )  (5) 

The  terms  and  are  given  by 

K2=//2max(v,+i.Vy,v/y_i) 

(0) 

f4=max(0,//4-K2) 

V  is  expressed  as 

^JPi+l-2pi+Pi-)\  (7) 

Pi+l  +  2pi+Pi-l 


The  unsteady  calculation  starts  from  steady  state 
solutions  from  a  circumferential  average  calculation. 
In  the  circumferential  average  calculation,  only  one 
blade  channel  per  blade  row  is  used.  The  calculated 
exit  condition  of  the  stator  form  the  inlet  unsteady 
boundary  condition  for  the  downstream  rotor  blade 
row.  On  the  other  way  around,  the  nonuniform  inlet 
flow  for  the  rotor  calculated  from  circumferential 
average  solution  forms  the  stator  unsteady  condition 
at  the  exit.  The  unsteady  interaction  between  the 
rotating  and  stationary  blade  rows  are  perturbed  with 
respect  to  the  mean  values  (the  steady  state  solutions). 
The  perturbed  flow  quantity  at  the  stator  exit  can  then 
be  expressed  as: 

(11) 


The  values  of  chosen  as  0.5  and  0.05 

respectively.  C  is  a  coefficient  that  has  an  impact  on 
the  stability  and  accuracy  of  the  solution.  To 
minimize  dissipation  within  viscous  boundary  layer, 
is  calculated  from : 


— (1+— ^ 
At^  At^j 


(8) 


where  (pj  is  the  stator  steady  state  solution  (obtained 

from  circumferential  average  calculation),  ^  is  the 
stator  pitchwise  mean  value.  The  flow  quantity  at  the 
rotor  inlet  would  be  the  steady  rotor  value  add  the 
perturbation  from  the  stator  exit  (equals  rotor 
pitchwise  mean  multiplied  by  stator  exit 
perturbation), 

)/*  "^  ( ’  '^s  (12) . 
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Similarily,  the  flow  quantity  at  the  stator  exit  can  be 
expressed  as  the  steady  stator  exit  value  plus  the 
pertubation  from  the  rotor  inlet, 

(13) 

Since  only  pitch  per  blade  row  is  used,  and  the  blade 
numbers  are  different,  the  perturbation  period  felt  by 
the  rotor  is  different  from  that  felt  by  the  stator.  Thus, 
one  must  keep  the  flow  data  of  the  neighboring  flow 
channels  in  memory.  Although  the  unsteady  flow 
quantity  calculated  from  eqs  (12)  and  (13)  are  not 
periodical  from  blade  to  blade,  and  also  are  varied 
from  time  to  time,  the  steady  state  solutions  for  the 
two  blade  rows  are  unchanged  within  one  blade 
passing  period.  It  must  be  emphasized  that  the  initial 
steady  state  solutions  are  obtained  from 
circumferential  average  calculations.  The  steady  state 
solutions  are  updated  by  a  time  average  method  after 
every  blade  passing  period  is  completed.  The  iteration 
process  is  continued  until  a  converged  solution  is 
achieved. 

3.  BOUNDARY  CONDITIONS 

The  boundary  conditions  for  turbomachinery  inlet, 
exit,  solid  surface,  periodic  boundaries  and  at  mixing 
planes  are  discussed  below. 

One-dimensional  Riemann  Invariant  non¬ 
reflection  boundary  condition  is  adopted  at  the  inlet 
boundary.  At  the  exit  boundary,  the  pressure  at  mid- 
span  is  prescribed.  Radial  equilibrium  equation  is 
used  to  calculate  the  pressures  along  the  exit  plane 
from  hub  to  tip.  At  the  solid  surfaces,  both  non-slip 
boundary  condition  and  non-penetration  condition  are 
considered.  That  is  to  force  the  relative  contravariant 
velocities  on  the  wall  surface  to  be  zero.  Periodic 
boundary  condition  is  imposed  on  the  blade-to-blade 
upstream  of  the  leading  edge  and  downstream  of  the 
trailing  edge  for  each  blade  row  to  ensure  a  correct 
converged  solution.  This  is  to  average  the 
extrapolated  lower  boundary  (at  j=l)  and  the 
extrapolated  upper  boundary  (at  j=jmax)  values. 

For  circumferential  average  calculation,  the 
treatment  at  the  mixing  plane  is  an  important  step. 
The  mixing  planes  are  located  approximately  midway 
between  two  adjacent  blade  rows.  The  mixing  planes 
for  circumferential  average  calculations  are  treated  as 
interior  points  rather  than  ordinary  inlet  or  exit 
boundary  points.  The  present  mixing  plane  treatment 


is  similar  to  that  of  Chen  [1996]  and  Denton  [1992]. 
The  variables  at  one  side  of  mixing  plane  are  first 
extrapolated  to  obtain  the  values  at  mixing  plane  and 
then  circumferentially  area  averaged. 

y=ymax  y=ymax 

<l>'ave=  S  («>/)4/(  E  4)  (14) 

j=\  J  J  j=\ 
y=ymax  j=jmsx 

/ave=  S  (<)4/(  S  4)  (15) 

j=\  J  ^  j^l 

The  superscripts  s  and  r  denote  stator  and  rotor 
respectively.  The  new  circumferential  average  values 
at  the  mixing  plane  are  the  simple  arithmetic  mean  of 
the  two  circumferentially  averaged  values  from  two 
sides. 

^ave  ~  ^^ave  ^ave^ ^ ^  (^^) 

The  variables  at  either  side  of  the  mixing  plane  are 
adjusted  proportionally  to  the  ratio  between  the  new 
and  old  circumferential  values. 

rave 

These  values  are  then  considered  as  the  new  boundary 
values  at  the  mixing  planes  for  the  next  step 
calculation.  This  mixing  plane  modeling  results  in 
circumferentially  nonuniform  flow  properties  at  the 
mixing  planes. 

Since  unequal  pitch  channels  are  used,  there  is  no 
unsteady  periodic  boundary  condition  as  that  in  a 
steady  state  sense.  The  unsteady  values  at  the  steady 
periodic  lines  are  simply  extrapolated  from  interior 
points  at  each  time  step. 

4.  RESULTS  AND  DISCUSSIONS 

The  UTRC  Large  Scale  Turbine  is  used  here  as  the 
baseline  configuration.  Its  geometrical  parameters 
are  shown  in  Table  1.  In  the  present  study,  only  the 
first  stator  and  rotor  are  considered.  The  gap  between 
stator  and  rotor  is  15%  stator  chord. 

Table  1  UTRC  11/2  Stage  Large  Scale  Turbine 
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1  st  stator 

1st  rotor 

2nd  stator 

blade  number 

22 

28 

22 

chord  (m) 

0.151 

0.161 

0.164 

hub  radius  (m) 

0.610 

0.610 

0.610 

tip  radius  (m) 

0.762 

0.762 

0.762 

The  grid  system  for  the  turbine  is  generated  using 
the  interactive  grid  generation  procedure  developed 
by  Lee  et  al.  [1995].  The  grid  system  is  the  same  as 
that  used  in  a  previous  study  [Chen  and  Lee  (1998)]. 
A  two-dimensional  view  of  the  the  grid  system  for  the 
1.5  stage  turbine  is  shown  in  fig.  1.  The  grid  system 
has  89x33x25  grids  for  the  first  stator,  85x33x25 
grids  for  the  rotor,  and  99x33x25  grids  for  the  second 
stator.  The  total  number  of  grids  is  225,225.  For  the 
present  study,  only  the  stator  and  rotor  grids  are 
considered.  Good  orthogonality,  smoothness  and 
clustemess  are  achieved. 

The  circumferential  average  solution  is  used  as  the 
initial  condition  for  unsteady  calculations.  Fig.  2 
shows  the  circumferential  average  solution  for  the 
stator  from  near  the  hub  to  near  near  the  tip.  The 
results  compared  well  with  the  data  [Dring  et  al. 
(1982)].  Fig.  3  shows  the  circumferential  average 
solution  for  the  rotor.  Some  uncertainty  in  the  data  for 
the  rotor  is  shown  in  the  figure.  Again  the  agreement 
is  good.  Fig.  4  is  the  pressure  distributions  at  various 
time  steps  for  the  stator  and  rotor  at  midspan.  The 
circles  in  the  figure  are  the  time  average  results  over 
one  blade  passing  period.  The  time  average  results 
become  the  new  steady  state  solution.  The  final 
converged  steady  state  solutions  at  several  blade 
heights  are  shown  in  figs.  2  and  3.  The  solutions  for 
both  the  initial  circumferential  average  solution  and 
the  final  time  averaged  steady  state  solutions  are 
similar  although  some  discrepancies  are  found  near 
the  hub  and  the  tip. 

The  pressures  at  6  locations  (3  near  stator  exit  and 
3  near  rotor  inlet  shown  in  fig.  5)  are  monitored  to 
make  sure  the  unsteady  pressures  do  reach  periodical 
solutions.  (Note  that  the  figure  is  plotted  of  equal 
pitch  used  in  the  other  calculations)  Fig.  6  shows  the 
unsteady  pressures  at  these  6  monitor  points.  It  can  be 
seen  that  the  unsteady  pressures  reach  near  periodical 
condition  very  quickly  after  only  one  blade  passing 
period. 

The  stator  exit  velocities  at  the  midspan  mixing 
plane  boundary  are  shown  in  fig.  7  for  one  rotor  blade 
passing  period.  The  stator  wake  at  pitch=:0.9  has  some 
significant  velocity  fluctuation  due  to  interaction  with 
the  downstream  rotor.  The  global  wave  moves  from 


pitch=l  to  pitch=0  due  to  rotor  motion.  The  time 
average  velocity  at  stator  exit  mixing  plane  is  shown 
in  fig.  8.  It  should  be  pointed  out  that  the  calculated 
result  has  been  shifted  tangentially  in  order  to 
compare  against  the  original  measurement  data. 
Although  the  predicted  wake  velocity  deficit  is  larger, 
the  overall  trend  agrees  quite  well  with  data.  Fig.  9  is 
the  rotor  inlet  velocities  at  the  mixing  plane 
boundary.  Global  perturbation  wave  motion  is 
moving  in  opposite  direction  relative  to  that  in  the 
stator.  The  stator  viscous  wake  influence  on  the  rotor 
is  obvious  as  it  moves  with  with  the  global  wave.  The 
results  show  that  the  rotor  influence  on  the  the 
upstream  stator  is  basically  a  potential  effect.  And, 
the  stator  influence  on  the  rotor  both  the  potential  and 
viscous  wake  effects  are  important.  The  time  average 
velocity  at  rotor  inlet  mixing  plane  is  shown  in  fig. 
10.  The  fluctuation  is  larger  than  that  due  to  the  stator 
wake  (fig.  8). 

The  pressure  fluctuation  amplitudes  for  the 
unsteady  calculations  are  shown  in  fig.  11.  The 
predictions  agreed  quite  well  with  the  data  especially 
the  maximum  amplitude  and  trend  at  the  critical  areas 
(trailing  edge  of  stator  and  leading  edge  of  rotor), 
although  some  overprediction  on  the  suction  surface 
of  stator  and  underprediction  on  the  suction  surface  of 
rotor. 

The  analysis  shows  that  a  three-dimensional 
unsteady  analysis  using  single  unequal  pitch  for  a 
turbine  stage  can  be  successful  if  the  interactions 
between  the  stationary  and  rotating  blade  rows  are 
properly  modeled.  The  biggest  advantage  of  using 
this  approach  is  computational  efficiency  and 
maintain  the  blade  configuration  intact.  Thus,  the 
predicted  flowfield  could  be  more  realistic.  The 
computations  were  performed  on  a  DEC-8400 
machine.  Roughly  1500  time  steps  per  stator  passing 
period,  and  3  subiterations  per  time  steps  were  used. 
One  entire  unsteady  calculation  (including  initial 
circumferential  average  calculation)  of  about  10  stator 
passing  periods  (or  roughly  13  rotor  passing  periods) 
took  about  6  days  in  CPU. 

5.  CONCLUSIONS 

A  loosely  coupled  blade  row  analysis  method  is 
used  on  a  three-dimensional  single  stage  turbine 
unsteady  analysis.  Unlike  the  other  loosely  coupled  or 
fully  coupled  methods,  the  present  approach  uses  only 
one  blade  pitch  for  each  blade  row  while  the  blade 
configuration  is  unmodified.  Circumferential  average 
analysis  result  using  the  same  computational  domain 
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is  used  as  the  initial  condition  for  unsteady  analysis. 
The  use  of  circumferential  average  solution  as  the 
starting  point  makes  the  predicted  unsteady  periodical 
solution  can  nearly  be  achieved  after  only  one  blade 
passing  period.  The  steady  state  solution  is  a  time 
average  result  over  a  complete  blade  passing  period. 
Basically  both  steady  state  solution  and 
circumferential  average  solution  agree  well  with 
measured  data,  although  some  differences  exist 
between  the  two  solutions.  The  predicted  unsteady 
solutions  show  thw  that  the  unsteady  interaction 
between  two  adjacent  blade  rows  can  be  reasonably 
well  simulated.  It  is  clear  that  the  influence  from  the 
downstream  blade  row  basically  is  due  to  a  potential 
effect  and  the  influence  from  the  upstream  blade  row 
is  due  both  to  viscous  wake  and  potential  effects. 
Despite  the  present  approach  is  considered  as  an 
approximate  solution  compared  with  a  fully  coupled 
unsteady  calculation,  the  result  is  very  encouraging. 
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Fig.  7  Unsteady  velocity  at  stator  exit. 
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Fig.  1 1  Unsteady  pressure  fluctuation  amplitudes. 
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